/*
APCED - Analysis for NH cohort
Hyesung "Hace" Oh
Started 11/07/2021
Last worked on: 03/14/2023
*/
****************************************************************************
*GENERAL*
clear all

*DATE*
global currdate: display %td_CCYY_NN_DD date(c(current_date), "DMY")
global date = subinstr(trim("$currdate"), " ", "_", .)

*LOG*
cd "P:\apced\shared\Aim2\source"
capture log close 
log using ..\output\analysis_NH_cohort_$date.log, text replace

*MAIN*
capture program drop main
program define main
	
	do ..\source\programs.do 
	use ..\source\final_nh_af.dta, clear
	capture drop _merge
	keep if bene_id_18900 != ""
	merge 1:1 bene_id_18900 using "P:\apced\shared\Data retrieval\adrd_ffs_cohort_newvars.dta", keepusing(er_count_last30 *_8mo_* *_30d_*) keep(3)
	alt_model_clean
	unadj_covar_assess
	gen_pinv_weights
	pinv_outcome_plots
	gen_pvol_weights
	areg_trt_models
	areg_alt_models //Run and store estimates 
	output_results_1
	var_scaling
	ipwra_models //Run and store estimates -- NEED TO WORK ON THESE
	output_results_2 // Output areg and ipwra models
	other_sensitivity //NEW AS OF JULY 1
	output_results_3
	
	log close

end



capture program drop alt_model_clean
program define alt_model_clean

	capture drop high_volume
	gen high_volume = (all_total >= 8)
	
	capture drop apc_volume
	gen apc_volume = .
	replace apc_volume = 0 if high_volume == 0 
	replace apc_volume = 1 if apc_rate < 0.5 & high_volume == 1
	replace apc_volume = 2 if apc_rate >= 0.5 & high_volume == 1
	
	label variable apc_volume "High volume and APC rate"
	capture label drop a_v
	label define a_v 0 "Low volume (< 8 E+M visits)" 1 ">= 8 E+M visits &<50% APC rate" 2 ">= 8 E+M visits &>= 50% APC rate"
	label values apc_volume a_v
	tab apc_volume
	
	capture drop state_group
	egen state_group = group(state)
	
	keep if apc_involve != . 
	
	*Generate an "all ER" variable
	capture drop all_er_count
		gen all_er_count = er_count + er_count_last30
		
	capture drop ever_base_er
		gen ever_base_er = er_count > 0
		tab ever_base_er
		
	capture drop er_base_decile
		xtile er_base_decile = all_er_count, nq(10)
		tab er_base_decile
		
	capture drop er_above_med
		sum er_count, detail
		local p75 = `r(p75)'
	capture drop er_above_p75
		gen er_above_p75 = (er_count >= `p75')
		tab er_above_p75
		
	capture drop tot_fac_phys_visits
		bys PROV0475: egen tot_fac_phys_visits = sum(phys_total)
	capture drop fac_any_md
		gen fac_any_md = tot_fac_phys_visits > 0
		
	capture drop fac_rurality_*
		tab fac_rurality, gen(fac_rurality_)
		
	capture drop hosptl_30day100
		gen hosptl_30day100 = hosptl_30day * 100
		
	capture drop hospice_death100
		gen hospice_death100 = hospice_death * 100
		
	capture drop imp_overall_rt
		gen imp_overall_rt = overall_rating_prev_impute
	

end

*NH quality covariates:  To control for the nursing home quality effects, you should include several nursing home characteristics such as staffing variables, star-ratings, profit status, chain status, size, share of Medicaid paid patients.

*facpoor - OSCAR: Low resource facility based on payer mix
*totbeds - OSCAR: Best-guess total beds in facility
*agg_adm - Number of Admissions
*paymcaid - OSCAR: Pct Medicaid primary payer
*multifac - OSCAR: Facility is part of a chain
*profit -  OSCAR: Facility is run for profit
*hospbase - OSCAR: Facility is hospital-based
*i.fac_rurality -  2009:Metro/Micro Indicator Code--0 = Not,1 = Metro,2 = Micro
*rn2nrs -  OSCAR: Ratio of RN/RN+LPN (0,1)
*rnhrppd - OSCAR: Total RN hrs/day/resident (adjusted)
*lpnhrppd - OSCAR: Total LPN hrs/day/resident (adjusted)
*cnahrppd - OSCAR: Total CNA hrs/day/resident (adjusted)

*facpoor totbeds agg_adm paymcaid multifac profit hospbase i.fac_rurality rn2nrs rnhrppd lpnhrppd cnahrppd anymdex

capture program drop unadj_covar_assess
program define unadj_covar_assess

	tab apc_involve


*agg_adm_prev_impute hospbase_prev_impute rnhrppd_prev_impute lpnhrppd_prev_impute cnahrppd_prev_impute overall_rating_prev_impute

capture drop agg_impute hospbase_impute rnhrppd_impute lpnhrppd_impute cnahrppd_impute overall_impute 
capture drop fac_micro

gen agg_impute = agg_adm_prev_impute
gen hospbase_impute = hospbase_prev_impute
gen rnhrppd_impute = rnhrppd_prev_impute
gen lpnhrppd_impute = lpnhrppd_prev_impute
gen cnahrppd_impute = cnahrppd_prev_impute
gen overall_impute = overall_rating_prev_impute
gen fac_micro = fac_metro != 1 & fac_not_metro_micro != 1

	foreach v in apc_involve apc_volume {

	*Assess covariates
	table1by_gen age adrd_years all_total female white black hispanic ccw dual_9mo er_count acute_days ip_days ltc_days hh_days hospc_days metro micro ami anem asth atrfib chf chrkid copd depr diab glauc hipfx hyprlip bph hyperten hypothy ischhd osteo arthr stroke schi schiot bipl pvd ulcers psds drug cancer facpoor totbeds paymcaid payother multifac profit hmo_share sdi_score agg_impute hospbase_impute rnhrppd_impute lpnhrppd_impute cnahrppd_impute overall_impute fac_metro fac_micro snf_bed hha hosptl_30day hospice_death using ..\output\table1_NH_by_`v', by(`v') mean_dec_place(2) sd_dec_place(2) replace

	}

end


capture program drop areg_trt_models
program define areg_trt_models

	tab high_volume apc_involve
	
	tab high_volume, sum(apc_rate)
	
	tab apc_involve, sum(high_volume)
	
	tab apc_involve
	
	areg apc_rate high_volume i.year i.state_group er_count acute_days age adrd_years female white black hispanic ccw dual_9mo ip_days ltc_days hh_days home_days hospc_days res_rurality ami anem asth atrfib chf chrkid copd depr diab glauc hipfx hyprlip bph hyperten hypothy ischhd osteo arthr stroke schi schiot bipl pvd ulcers psds drug cancer facpoor totbeds paymcaid payother multifac profit hmo_share sdi_score agg_adm_prev_impute hospbase_prev_impute rnhrppd_prev_impute lpnhrppd_prev_impute cnahrppd_prev_impute overall_rating_prev_impute i.fac_rurality snf_bed hha, absorb(hkzip) cluster(hkzip)

	

end


capture program drop areg_alt_models
program define areg_alt_models

tab apc_involve

eststo clear

display "APC involve as treatment; High volume as control"
display "Then, high volume as treatment"
foreach var of varlist hosptl_30day100 hospice_death100 {
	
		eststo, prefix(`var'_areg_inv_): areg `var' i.apc_involve i.year i.state_group er_count acute_days age adrd_years female white black hispanic ccw dual_9mo ip_days ltc_days hh_days home_days hospc_days i.res_rurality ami anem asth atrfib chf chrkid copd depr diab glauc hipfx hyprlip bph hyperten hypothy ischhd osteo arthr stroke schi schiot bipl pvd ulcers psds drug cancer facpoor totbeds paymcaid payother multifac profit hmo_share sdi_score agg_adm_prev_impute hospbase_prev_impute rnhrppd_prev_impute lpnhrppd_prev_impute cnahrppd_prev_impute overall_rating_prev_impute i.fac_rurality snf_bed hha, absorb(hkzip) cluster(hkzip)
		
		
		eststo, prefix(`var'_areg_vol_): areg `var' i.apc_volume i.year i.state_group er_count acute_days age adrd_years female white black hispanic ccw dual_9mo ip_days ltc_days hh_days home_days hospc_days i.res_rurality ami anem asth atrfib chf chrkid copd depr diab glauc hipfx hyprlip bph hyperten hypothy ischhd osteo arthr stroke schi schiot bipl pvd ulcers psds drug cancer facpoor totbeds paymcaid payother multifac profit hmo_share sdi_score agg_adm_prev_impute hospbase_prev_impute rnhrppd_prev_impute lpnhrppd_prev_impute cnahrppd_prev_impute overall_rating_prev_impute i.fac_rurality snf_bed hha, absorb(hkzip) cluster(hkzip)
		

}
	
end


capture program drop output_results_1
program define output_results_1

	foreach var of varlist hosptl_30day100 hospice_death100 {
		
		esttab `var'_areg_* using "../output/`var'_all_areg_NH.rtf", replace drop(*state_group* *year* *fac_rurality* er_count acute_days age adrd_years female white black hispanic ccw dual_9mo ip_days ltc_days hh_days home_days hospc_days *res_rurality* ami anem asth atrfib chf chrkid copd depr diab glauc hipfx hyprlip bph hyperten hypothy ischhd osteo arthr stroke schi schiot bipl pvd ulcers psds drug cancer facpoor totbeds paymcaid payother multifac profit hmo_share sdi_score agg_adm_prev_impute hospbase_prev_impute rnhrppd_prev_impute lpnhrppd_prev_impute cnahrppd_prev_impute overall_rating_prev_impute snf_bed hha) b(%9.1fc) ci(%9.1fc) par scalars(F rmse) sfmt(%9.6fc) mtitles("APC Involve FS" "APC/Volume-Int FS")  label nogaps compress 
		
	}

end


capture program drop var_scaling
program define var_scaling


foreach v in er_count acute_days age adrd_years female white black hispanic ccw dual_9mo ip_days ltc_days hh_days home_days hospc_days ami anem asth atrfib chf chrkid copd depr diab glauc hipfx hyprlip bph hyperten hypothy ischhd osteo arthr stroke schi schiot bipl pvd ulcers psds drug cancer facpoor totbeds paymcaid payother multifac profit hmo_share sdi_score agg_adm_prev_impute hospbase_prev_impute rnhrppd_prev_impute lpnhrppd_prev_impute cnahrppd_prev_impute overall_rating_prev_impute  snf_bed hha {
    
	quietly sum `v'
	if `r(max)'>0 & `r(max)' !=. {
		replace `v' = `v'/`r(max)'
	}
	
	
}


end

capture program drop ipwra_models
program define ipwra_models

eststo clear
	
foreach var of varlist hosptl_cnt er_cnt {
	eststo, prefix(`var'_ipwra_): teffects ipwra (`var' i.year i.state_group er_count acute_days age adrd_years female white black hispanic ccw dual_9mo ip_days ltc_days hh_days home_days hospc_days i.res_rurality ami anem asth atrfib chf chrkid copd depr diab glauc hipfx hyprlip bph hyperten hypothy ischhd osteo arthr stroke schi schiot bipl pvd ulcers psds drug cancer facpoor totbeds paymcaid payother multifac profit hmo_share sdi_score agg_adm_prev_impute hospbase_prev_impute rnhrppd_prev_impute lpnhrppd_prev_impute cnahrppd_prev_impute overall_rating_prev_impute i.fac_rurality snf_bed hha) (apc_involve i.year i.state_group er_count acute_days age adrd_years female white black hispanic ccw dual_9mo ip_days ltc_days hh_days home_days hospc_days i.res_rurality ami anem asth atrfib chf chrkid copd depr diab glauc hipfx hyprlip bph hyperten hypothy ischhd osteo arthr stroke schi schiot bipl pvd ulcers psds drug cancer facpoor totbeds paymcaid payother multifac profit hmo_share sdi_score agg_adm_prev_impute hospbase_prev_impute rnhrppd_prev_impute lpnhrppd_prev_impute cnahrppd_prev_impute overall_rating_prev_impute i.fac_rurality snf_bed hha)
	
		eststo, prefix(`var'_ipwra_): teffects ipwra (`var' i.year i.state_group er_count acute_days age adrd_years female white black hispanic ccw dual_9mo ip_days ltc_days hh_days home_days hospc_days i.res_rurality ami anem asth atrfib chf chrkid copd depr diab glauc hipfx hyprlip bph hyperten hypothy ischhd osteo arthr stroke schi schiot bipl pvd ulcers psds drug cancer facpoor totbeds paymcaid payother multifac profit hmo_share sdi_score agg_adm_prev_impute hospbase_prev_impute rnhrppd_prev_impute lpnhrppd_prev_impute cnahrppd_prev_impute overall_rating_prev_impute i.fac_rurality snf_bed hha) (apc_volume i.year i.state_group er_count acute_days age adrd_years female white black hispanic ccw dual_9mo ip_days ltc_days hh_days home_days hospc_days i.res_rurality ami anem asth atrfib chf chrkid copd depr diab glauc hipfx hyprlip bph hyperten hypothy ischhd osteo arthr stroke schi schiot bipl pvd ulcers psds drug cancer facpoor totbeds paymcaid payother multifac profit hmo_share sdi_score agg_adm_prev_impute hospbase_prev_impute rnhrppd_prev_impute lpnhrppd_prev_impute cnahrppd_prev_impute overall_rating_prev_impute i.fac_rurality snf_bed hha)
		


}


end


capture program drop output_results_2
program define output_results_2

	foreach var of varlist hosptl_cnt er_cnt {

		
		esttab `var'_ipwra_* using "../output/`var'_ipwra_nh_models.rtf", replace drop(*year* *state_group* *fac_rurality* er_count acute_days age adrd_years female white black hispanic ccw dual_9mo ip_days ltc_days hh_days home_days hospc_days *res_rurality* ami anem asth atrfib chf chrkid copd depr diab glauc hipfx hyprlip bph hyperten hypothy ischhd osteo arthr stroke schi schiot bipl pvd ulcers psds drug cancer facpoor totbeds paymcaid payother multifac profit hmo_share sdi_score agg_adm_prev_impute hospbase_prev_impute rnhrppd_prev_impute lpnhrppd_prev_impute cnahrppd_prev_impute overall_rating_prev_impute snf_bed hha) b(%9.1fc) ci(%9.1fc) par scalars(F rmse) sfmt(%9.6fc) mtitles("`var': IPW RA w/ APC Involvement Trt" "`var': IPW RA w/ APC/Volume Trt") label nogaps compress
		
	}

end

capture program drop other_sensitivity
program define other_sensitivity 

	eststo clear

display "APC involve as treatment"
display "Then, APC/Volume as treatment"
display "Another alternative specification: High APC x High Volume interaction"

	capture drop apc_rate50
	gen apc_rate50 = apc_rate >= 0.5
	
	capture drop low_volume 
	gen low_volume = high_volume == 0
	
	capture drop highvol_apc50 
	gen highvol_apc50 = .
	replace highvol_apc50 = 0 if apc_rate50 == 0 & high_volume == 0
	replace highvol_apc50 = 1 if apc_rate50 == 1 & high_volume == 0
	replace highvol_apc50 = 2 if high_volume == 1
	
	capture drop lowvol_apc50_int
	gen lowvol_apc50_int =  low_volume * apc_rate50

	capture drop violator_*
	
		foreach var of varlist hosptl_30day100 hospice_death100 {

		eststo, prefix(`var'_areg_int): areg `var' i.highvol_apc50 i.year i.state_group er_count acute_days age adrd_years female white black hispanic ccw dual_9mo ip_days ltc_days hh_days home_days hospc_days i.res_rurality ami anem asth atrfib chf chrkid copd depr diab glauc hipfx hyprlip bph hyperten hypothy ischhd osteo arthr stroke schi schiot bipl pvd ulcers psds drug cancer facpoor totbeds paymcaid payother multifac profit hmo_share sdi_score agg_adm_prev_impute hospbase_prev_impute rnhrppd_prev_impute lpnhrppd_prev_impute cnahrppd_prev_impute overall_rating_prev_impute i.fac_rurality snf_bed hha, absorb(hkzip) cluster(hkzip)
		
		eststo, prefix(`var'_ipw_int): teffects ipwra (`var' i.year i.state_group er_count acute_days age adrd_years female white black hispanic ccw dual_9mo ip_days ltc_days hh_days home_days hospc_days i.res_rurality ami anem asth atrfib chf chrkid copd depr diab glauc hipfx hyprlip bph hyperten hypothy ischhd osteo arthr stroke schi schiot bipl pvd ulcers psds drug cancer facpoor totbeds paymcaid payother multifac profit hmo_share sdi_score agg_adm_prev_impute hospbase_prev_impute rnhrppd_prev_impute lpnhrppd_prev_impute cnahrppd_prev_impute overall_rating_prev_impute i.fac_rurality snf_bed hha) (highvol_apc50 i.year i.state_group er_count acute_days age adrd_years female white black hispanic ccw dual_9mo ip_days ltc_days hh_days home_days hospc_days i.res_rurality ami anem asth atrfib chf chrkid copd depr diab glauc hipfx hyprlip bph hyperten hypothy ischhd osteo arthr stroke schi schiot bipl pvd ulcers psds drug cancer facpoor totbeds paymcaid payother multifac profit hmo_share sdi_score agg_adm_prev_impute hospbase_prev_impute rnhrppd_prev_impute lpnhrppd_prev_impute cnahrppd_prev_impute overall_rating_prev_impute i.fac_rurality snf_bed hha), osample(violator_vol_`var')
	
	}

end

capture program drop output_results_3
program define output_results_3

	foreach var of varlist hosptl_30day100 hospice_death100 {
		
		esttab `var'_areg_int* using "../output/`var'_NH_areg_int_models.rtf", replace drop(*state_group* *year* er_count acute_days age adrd_years female white black hispanic ccw dual_9mo ip_days ltc_days hh_days home_days hospc_days *res_rurality* ami anem asth atrfib chf chrkid copd depr diab glauc hipfx hyprlip bph hyperten hypothy ischhd osteo arthr stroke schi schiot bipl pvd ulcers psds drug cancer facpoor totbeds paymcaid payother multifac profit hmo_share sdi_score agg_adm_prev_impute hospbase_prev_impute rnhrppd_prev_impute lpnhrppd_prev_impute cnahrppd_prev_impute overall_rating_prev_impute *fac_rurality* snf_bed hha) b(%9.1fc) ci(%9.1fc) par scalars(F rmse) sfmt(%9.6fc) mtitles("`var': AReg Int Trt") label nogaps compress
		
		esttab `var'_ipw_int* using "../output/`var'_NH_ipw_int_models.rtf", replace drop(*state_group* *year* er_count acute_days age adrd_years female white black hispanic ccw dual_9mo ip_days ltc_days hh_days home_days hospc_days *res_rurality* ami anem asth atrfib chf chrkid copd depr diab glauc hipfx hyprlip bph hyperten hypothy ischhd osteo arthr stroke schi schiot bipl pvd ulcers psds drug cancer facpoor totbeds paymcaid payother multifac profit hmo_share sdi_score agg_adm_prev_impute hospbase_prev_impute rnhrppd_prev_impute lpnhrppd_prev_impute cnahrppd_prev_impute overall_rating_prev_impute *fac_rurality* snf_bed hha) b(%9.1fc) ci(%9.1fc) par scalars(F rmse) sfmt(%9.6fc) mtitles("`var': IPWRA Int Trt") label nogaps compress
		
	}

end

capture program drop gen_pinv_weights
program define gen_pinv_weights


	capture drop pinv1
	capture drop pinv2
	capture drop pinv3
	mlogit apc_involve i.year i.state_group er_count acute_days age adrd_years female white black hispanic ccw dual_9mo ip_days ltc_days hh_days home_days hospc_days i.res_rurality ami anem asth atrfib chf chrkid copd depr diab glauc hipfx hyprlip bph hyperten hypothy ischhd osteo arthr stroke schi schiot bipl pvd ulcers psds drug cancer facpoor totbeds paymcaid payother multifac profit hmo_share sdi_score agg_adm_prev_impute hospbase_prev_impute rnhrppd_prev_impute lpnhrppd_prev_impute cnahrppd_prev_impute overall_rating_prev_impute i.fac_rurality snf_bed hha, base(0) robust difficult	
	predict pinv*
	sum pinv1 pinv2 pinv3
	
	
graph box pinv1, over(apc_involve, label(labsize(small))) noout	///
title("{bf:NH cohort}", color(black) size(medsmall))	///
subtitle("Multinomial regression: bene/NH chars, state FEs", size(small)) title("{bf:Figure A7.} Predicted probabilities of being in the minimal APC involvement group")	///
ytitle("Predicted probability of being in" "{bf:minimal} APC treatment group", size(medsmall)) ysize(9) xsize(16)	///
ylabel(#10, angle(0) labsize(small))	///
ysc(titlegap(2))	///
box(1, color(dkgreen))	///
graphregion(color(white))	///
name(nhinv_p1, replace)
	graph export ..\output\FigA7.tif, replace width(4000)
	
graph box pinv2, over(apc_involve, label(labsize(small))) noout	///
title("{bf:NH cohort}", color(black) size(medsmall)) title("{bf:Figure A8.} Predicted probabilities of being in the moderate APC involvement group")	///
subtitle("Multinomial regression: bene/NH chars, state FEs", size(small))	///
ytitle("Predicted probability of being in" "{bf:moderate} APC treatment group", size(medsmall)) ysize(9) xsize(16)	///
ylabel(#10, angle(0) labsize(small))	///
ysc(titlegap(2))	///
box(1, color(dkgreen))	///
graphregion(color(white))	///
name(nhinv_p2, replace)
	graph export ..\output\FigA8.tif, replace width(4000)


graph box pinv3, over(apc_involve, label(labsize(small))) noout	///
title("{bf:NH cohort}", color(black) size(medsmall))	///
subtitle("Multinomial regression: bene/NH chars, state FEs", size(small)) title("{bf:Figure A9.} Predicted probabilities of being in the extensive APC involvement group")	///
ytitle("Predicted probability of being in" "{bf:extensive} APC treatment group", size(medsmall)) ysize(9) xsize(16)	///
ylabel(#10, angle(0) labsize(small))	///
ysc(titlegap(2))	///
box(1, color(dkgreen))	///
graphregion(color(white))	///
name(nhinv_p3, replace)
	graph export ..\output\FigA9.tif, replace width(4000)


end


capture program drop pinv_outcome_plots
program define pinv_outcome_plots

	capture drop pinv_tails
	sum pinv1, detail
	gen pinv_tails = (pinv1 < r(p5) | pinv1 > r(p95))
	
	
	graph twoway lpoly hosptl_30day pinv1 if apc_involve==0 & pinv_tails == 0 || lpoly hosptl_30day pinv1 if apc_involve==1 & pinv_tails == 0 || lpoly hosptl_30day pinv1 if apc_involve==2 & pinv_tails == 0, scheme(sj) legend(pos(1) ring(0) col(1) lab(1 "Low APC") lab(2 "Limited APC") lab(3 "Extensive APC")) xtitle("Predicted probability: Low APC") ytitle("Probability: Hospitalized during final 30 days") title("Pr(30-day hospitaliztion) vs. Predicted Pr(Low APC)")
	graph export ..\output\NH_Hosp30_PropLowAPC.png, replace

	
	graph twoway lpoly hospice_death pinv1 if apc_involve==0 & pinv_tails == 0 || lpoly hospice_death pinv1 if apc_involve==1 & pinv_tails == 0 || lpoly hospice_death pinv1 if apc_involve==2 & pinv_tails == 0, scheme(sj) legend(pos(1) ring(0) col(1) lab(1 "Low APC") lab(2 "Limited APC") lab(3 "Extensive APC")) xtitle("Predicted probability: Low APC") ytitle("Probability: Death while enrolled in hospice") title("Pr(Death enrolled in hospice) vs. Predicted Pr(Low APC)")
	graph export ..\output\NH_HospiceDeath_PropLowAPC.png, replace

	
end


capture program drop gen_pvol_weights
program define gen_pvol_weights


	capture drop pvol1
	capture drop pvol2
	capture drop pvol3
	mlogit apc_volume i.year i.state_group er_count acute_days age adrd_years female white black hispanic ccw dual_9mo ip_days ltc_days hh_days home_days hospc_days i.res_rurality ami anem asth atrfib chf chrkid copd depr diab glauc hipfx hyprlip bph hyperten hypothy ischhd osteo arthr stroke schi schiot bipl pvd ulcers psds drug cancer facpoor totbeds paymcaid payother multifac profit hmo_share sdi_score agg_adm_prev_impute hospbase_prev_impute rnhrppd_prev_impute lpnhrppd_prev_impute cnahrppd_prev_impute overall_rating_prev_impute i.fac_rurality snf_bed hha, base(0) robust difficult	
	predict pvol*
	sum pvol1 pvol2 pvol3
	
	
graph box pvol1, over(apc_volume, label(labsize(vsmall))) noout	///
title("{bf:NH cohort}", color(black) size(medsmall))	///
subtitle("Multinomial regression: bene/NH chars, state FEs", size(small)) title("{bf:Figure A10.} Predicted probabilities of being in the low volume group")	///
ytitle("Predicted probability of being in" "{bf:low} APC vol group", size(medsmall)) ysize(9) xsize(16)	///
ylabel(#10, angle(0) labsize(small))	///
ysc(titlegap(2))	///
box(1, color(dkgreen))	///
graphregion(color(white))	///
name(nhvol_p1, replace)
	graph export ..\output\FigA10.tif, replace width(4000)

	
graph box pvol2, over(apc_volume, label(labsize(vsmall))) noout	///
title("{bf:NH cohort}", color(black) size(medsmall))	///
subtitle("Multinomial regression: bene/NH chars, state FEs", size(small)) title("{bf:Figure A11.} Predicted probabilities of being in the high volume + <50% APC involvement group")	///
ytitle("Predicted probability of being in" "{bf:high} E+M <50% APC group", size(medsmall)) ysize(9) xsize(16)	///
ylabel(#10, angle(0) labsize(small))	///
ysc(titlegap(2))	///
box(1, color(dkgreen))	///
graphregion(color(white))	///
name(nhvol_p2, replace)
	graph export ..\output\FigA11.tif, replace width(4000)
	


graph box pvol3, over(apc_volume, label(labsize(vsmall))) noout	///
title("{bf:NH cohort}", color(black) size(medsmall))	///
subtitle("Multinomial regression: bene/NH chars, state FEs", size(small)) title("{bf:Figure A12.} Predicted probabilities of being in the high volume + >50% APC involvement group")	///
ytitle("Predicted probability of being in" "{bf:high} E+M 50%+ APC group", size(medsmall))	ysize(9) xsize(16) ///
ylabel(#10, angle(0) labsize(small))	///
ysc(titlegap(2))	///
box(1, color(dkgreen))	///
graphregion(color(white))	///
name(nhvol_p3, replace)
	graph export ..\output\FigA12.tif, replace width(4000)

	

end





main





